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ABSTRACT 

We provide computationally convenient expressions for all marginal distributions of the polarization 
CMB power spectrum distribution P{Ci\(Ji), where Ci = {Cj'^ , Cj^ , Cf^ , Cf^} denotes the set of 
ensemble averaged polarization CMB power spectra, and at = {trj"^, u™, cr™, ^} the set of the 
realization specific polarization CMB power spectra. This distribution describes the CMB power 
spectrum posterior for cosmic variance limited data. The expressions derived here are general, and may 
be useful in a wide range of applications. Two specific applications are described in this paper. First, 
we employ the derived distributions within the CMB Gibbs sampling framework, and demonstrate a 
new conditional CMB power spectrum sampling algorithm that allows for different binning schemes 
for each power spectrum. This is useful because most CMB experiments have very different signal-to- 
noise ratios for temperature and polarization. Second, we provide new Blackwell-Rao estimators for 
each of the marginal polarization distributions, which are relevant to power spectrum and likelihood 
estimation. Because these estimators represent marginals, they are not affected by the exponential 
behaviour of the corresponding joint expression, but converge quickly. 

Subject headings: cosmic microwave background — cosmology: observations — methods: numerical 



1. INTRODUCTION 

During the last few decades cosmology has evolved 
from a data starved branch of astrophysics, into a data 
driven high-precision science in which theories may be 
subjected to stringent observational tests. This revolu- 
tion has to a large extent been driven by steadily improv- 
ing observations of the cosmic microwave background 
(CMB), allowing cosmologists to have a close-up look 
at the very young universe . Two leading experiments 

and WMAP 



were the COBE -DMR (Sm oot et al.l 11992 
([Bennett et al.lf2003f ) satellite missions, while the third 
generation experiment, Planck, will be launched late this 
year. 

As observations continue to improve, increasingly de- 
manding requirements are imposed on the data analy- 
sis. While rather crude approximations may be accept- 
able when interpreting low signal-to-noise data, the situ- 
ation is very different in the mid and high signal-to-noise 
regime. Here, even "small" effects become clearly visible, 
and may potentially compromise any cosmological con- 
clusion. Using accurate methods in this regime is criti- 
cal. Some real- world issues relevant to the CMB problem 
are non-cosmological foregrounds, improper noise and/or 
beam characterization, and sub-optimal likelihood ap- 
proximations. 

In 2004, a new app roach to CMB analysis was p roposed 
and im plem ented by I Jewell et a .1 ([2004) . Wandelt et al.l 
(|2004f ) and lEriksen et all (|2004( l. Rather than taking 
the traditional ap proximate Monte Carlo approach (e.g., 
iHivon et al.l 12003 ). this new method employs the Gibbs 
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sampling algorithm to facilitate exact (in the maximum- 
likelihood sense), global and efficient analysis of even 
high-resolution data sets. Equally important, the Gibbs 
sampling framework has unique capabilities for error 
propagation, as it allows for easy marginalization over 
virtually any auxiliary stochastic field. One important 
example is that of non-cosmological foregrounds. 

Since then, the method has been generaliz ed to han- 
dle polarized CMB data (ILarson et al.l 120 07^ and joint 
foreground and CMB analvsis (|Eriksen et a l. 2008a), and 
has been applied mos t successfully to the W MAP data 
(lO'Dwver et all I200I lEriksen et al"] l2007al lbl. [2008b). 
Some useful examples of issues correctly identified by the 
Gibbs sampler, but missed by other techni ques, are 1) the 
first-y e ar WMAP likelihood bias at £ < 30 ([Eriksen et~all 
l2007at iHinshaw et al.ll2007f ). 2) foreground resid uals in 
the 3-y ear WMAP polarization sky maps (Erikse n et al.l 
l2007 m. and 3) residual monopole and dip ole compo- 
nents in the 3-year ternperat ure sky maps (jEriksen et al.l 
I2008bt IHinshaw etall I2OO8O . Following up on these 
methodological advances, the WMAP team adopted the 
Gibbs sampler as a central component in their analysis 
of the 5-year data, and, in fact, their default low-£ tem- 
perature likelihood module is precisely the G ibbs-based 
Blackw ell-Rao code written and published by lChu et al.l 
(|2005h . 

While WMAP has done an excellent job on character- 
izing the large-scale CMB temperature fluctuations, the 
current frontier in CMB science is polarization. In just 
a few years, full-sky high-sensitivity data will be avail- 
able from Planck. And then, very likely, the situation 
will be quite analogous to the one WMAP experienced 
in the temperature case: Having robust, exact methods 
that allows for proper characterization and propagation 
of systematics will be essential in the mid to high signal- 
to-noise regime. The Gibbs sampler is among the leading 
candidates to serve such a purpose. 

Unfortunately, the Gibbs sampler, as currently de- 
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scribed in the literature, has two major hmitations that 
needs to be resolved before this promise can be fulfilled. 
First, the direct Gibbs sampler is inherently inefRcient in 
the low signal-to-noise regime, because the step size be- 
tween two consecutive samples is determined by cosmic 
variance alone, whereas the full posterior width is deter- 
mined by noise. Second, it is non-trivial to establish a 
full likelihood approximation from the samples produced 
by the Gibbs sampler, because of the dimensionality of 
the underlying space. Both of these issues are currently 
under development, and reports are expected in the near 
future (Jewell et al. 2008; Rudjord et al. and Eriksen et 
al., in preparation). 

In the present paper, we take a small but important 
first step towards resolving these issues, by considering 
the marginal and conditional densities of the probabil- 
ity distribution P{Cg\ai), where is the ensemble aver- 
aged CMB power spectrum, and ai is the observed power 
spectrum of one given CMB realization. This distribu- 
tion plays a crucial role within the CMB Gibbs sampling 
framework. On the one hand, it forms one of the two con- 
ditionals in the main sampling scheme. On the other, it 
is the kernel of the Blackwell-Rao estimator. Being able 
to describe this analytically in different forms is therefore 
very useful. Two specific applications are demonstrated 
in this paper, namely 1) a Ci sampling algorithm that 
allows for different binning schemes in each of the po- 
larization components, and 2) Blackwell-Rao estimators 
for each oi P{Cj^\Cj^ ,Cf^ ,d), P{Cj^\Cf^,d), and 
P{Cf^\d). Further applications will be demonstrated in 
the papers mentioned above. 

We also note that these expressions are completely gen- 
eral, and may prove useful for any other method that 
considers both Ce and the CMB sky signal s as free vari- 
ables. One such example i s the CMB Hamilto nian sam- 
pler recently developed bv lTavlor et al.l ()2007l ). 

2. NOTATION AND DATA MODEL 

We now introduce a statistical model for the CMB ob- 
servations, and define our notation. First, we assume 
that the data may be modelled by a signal and a noise 
term. 



d(n) = s{n) + n{n). 



(1) 



Here, d is a 3-component (T, Q, U) Stokes' parameter 
vector observed in direction fi. s and n denote similar 
vectors, describing the CMB field and instrumental noise, 
respectively. Both the signal and the noise are assumed 
to be Gaussian distributed with zero mean and covari- 
ances C and N, respectively. (Note that we, for nota- 
tional simplicity, neglect real-world complications such as 
instrumental beams, frequency dependent observations 
or foreground components in this expression; the topic 
of this paper is the probability distribution P{Ci, s), and 
for this, all such issues are irrelevant.) 

Next, we additionally assume the CMB field to be sta- 
tistically isotropic. It is therefore useful to decompose 
the (T, Q, U) field into spin-weighted spherical harmon- 
ics (see, e.g., Zaldarriaga and Seljak 1997 for full details). 



with coefficients (a 



Because the spherical harmonics are orthogonal on the 
full sky, and B has opposite parity of T and E, the har- 



monic space CMB covariance matrix is given by 




(2) 
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In this expression, brackets denote ensemble averages, 
and the power spectrum therefore corresponds to 
a theory spectrum, similar to that produced by a 
Boltz mann code such as CMBFast (jSeliak fc Zaldarriagal 
Il996l ) . Note that Ct denotes the matrix of all power spec- 
tra, while a single component is indicated by subscripts 

(e.g.,cr). 

One can also define the realization specific power spec- 
trum, (T£, which is simply the averaged power in each 
multipolc for one given realization. 
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Explicitly, Ci is the power spectrum corresponding to 
some cosmological model, and ai is the power spectrum 
of one realization drawn from that model. It may there- 
fore be useful to imagine that observations of the CMB 
sky provide us with cr^, and from this we seek to con- 
strain the underlying cosmological theory, parametrized 
by Ci and summarized by the conditional distribution 

P{CAai). 

With this notation, it is straightforward to write down 
the joint probability distribution for the CMB sky signal, 
s, the CMB power spectrum, Ci and the data, d, 

P(s,Q,d) = P(d|s,Q)P(s,Q) 

^g-i(d-srN-(d-s)p(g^C',). (4) 

and the CMB posterior distribution, 
F(s,Q,d) 



P(s,Q|d) 



P(d) 



«P(s,C,,d) 



(5) 



These expressions involve two factors, namely the = 
(d - s)'^N"i(d - s) and the CMB probability distribu- 
tion P{s,Ce). Since we assume that the CMB field is 
isotropic and Gaussian, as discussed above, the latter 
may be written as (e.g., Larson et al. 2007) 

Pis,Ce)^Pis\Ci)P{Ce) 




n 



\Ci\' 



P{Ci) (6) 



\{P{a,\Ct)P{Ci), 



where P{Ci) is a prior on C^. P{ai\Ci) is recognized 
as an inverse Wishart distribution when interpreted as a 
function of C^. 

Before turning to the main topic of this paper, we 
recall that, for a probability distribution P(x,y), the 
marginal distribution is defined by P{x) = J P{x,y)dy, 
and the conditional by P{y\x) = P{x,y)/P{x). From 
these, one may also trivially derive Bayes' theorem, 
P{x\y) =^ P{y\x)P{x) / P[y). Therefore, for uniform pri- 
ors on both Cf, and s, which we assume in this paper, 
P{ai\Ci) cx P{Ci\oi). 

3. CMB POWER SPECTRUM DISTRIBUTIONS 

The main goal of this paper is to derive explicit expres- 
sions for the marginals, and thereby the conditionals, of 
P{Ct\crt). These are summarized in Table[T] 

As seen above, P{Ce,\ai) is given by the inverse Wishart 
distribution, which, in n dimensions and inclu ding the 
full normalization factor (jGupta &: Nagaijr200ClD . reads 



P{CA<yd 



(7) 



r„(^)|Q| 

Here r„ is the multivariate Gamma function. 

However, as discussed in Section [21 wc are in this paper 
interested in the special case for which Cf^ — Cf^ — 
0. In this case, the trace in equation [7] expands into a 
sum of two terms, and the Ci determinant factorizes into 
the product of a two-dimensional (T, E) determinant and 
Cf^. Thus, the joint distribution factorizes as 

P{CA<yi) - P{Cr. Cf^\a,)P{Cf^\<j,). (8) 

That is, Cf^ is independent of {Cj'^ ,Cj^ .Cf'^), and 
follows a one-dimensional inverse Wishart (or inverse 
Gamma) distribution. We will therefore not consider 
the BB component further in this paper. However, we 
note that if one is interested in exotic models for which 
{TB, EB} 7^ 0, the expressions derived in this paper will 
have to be revised accordingly. 

3.1. The {TE,EE) distribution, P{Cj^ ,Cf'^\ai) 

We start by considering the two-dimensional marginal 
distribution P{Cj^ , Cf^\(Jt), which is obtained by inte- 
grating PiCj^ , Cj^ , Cf^\cT,) over Cj^ , 



p{cr,crw^)^ / p{cr,cr,crwddc^, 

2^+1 



'TT 



oc \ai\ 
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(9) 



2 (^TT[^EE_((7TE\2 T^TT 

Note that the lower limit in this integral is defined by 
\Ce\ > 0, since the power spectrum covariance matrix in 
equation [2] must be positive definite. 
If we now define 



EE (C™)^ 
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and make the change of variable 

(2£+ l)/c 



y 



(10) 



(11) 



this expression is transformed into 
I 

p{cr,cr\a,)^^g^ 



^EE 
2 TTETE" 



(12) 

The integral in this expression is simply the Gamma 
function, 



2£- 1 
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y " e' 



'dy, 



(13) 



and, for our purposes, an irrelevant numerical normal- 
ization factor. Thus, the final distribution reads 
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(14) 



EE(C'TE-;2 



2aTEcTE + arCf^ 



Note that because TT and EE occur symmetri- 
cally in equation [7l the corresponding expression for 
P{Cj'^ , Cj^\ae) is obtained simply by interchanging EE 
and TT in equation [Til 



3.2. The {TT, EE) distribution, P{Cj ^ , C^ 



lEE 



Next, we consider P[Cj^ ,Cf^\ai), which is obtained 
by integrating P{Cj'^ ,Cj^ ,Cf^\ae) over Cj^ . Unfor- 
tunately, this distribution does not have a closed expres- 
sion, but can instead be written on the form 



P{Cl\Cr\a,)^ 



Here Ii denotes the integral 
»1 

h - 



h{l,A,B). (15) 



-1 (1 — x-^) 2 



■dx. 



(16) 



and we have defined the two dimensionless auxiliary pa- 
rameters 



A = 



B 
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(17) 
(18) 



However, the fact that I\ depends only on two dimen- 
sionless parameters and I, implies that it can easily be 
tabulated (and optionally splined for higher accuracy) 
for each I, and thus computationally efficient lookup- 
tables may be constructed. In most practical appli- 
cations, which typically require repeated evaluations of 
P{Cj'^ ,Cf^\at), Equation [m is therefore as useful for 
(Cj^, Cf^) as EquationHHis for {Cj^ ,Cf^), although 
implementationally a little more complicated. 

3.3. The marginal EE distribution, P{Cf^\ai) 

We now compute the corresponding one-dimensional 
marginals, and begin with P{Cf^\ae), by integrating 
equation [Til over C™. This is simplified by introduc- 
ing the new variable 



y 



„EE/^TE 



^TE/^EE 



'<EE 



(19) 
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TABLE 1 

CMB POWER SPECTRUM DISTRIBUTIONS 



Distribution Expression 



Joint {T, E} distribution 



P{Cj'r,Cj^,Cf^\ai) 


oc 


2l+l « 


Biveiriate marginals 



P{Cr,Cj^\ai) 
P{Cr,Cf^\ai) 



P{,Cr,Cf^\ai) 



We\' 



2£+l 



it-i 6 



(aTT(cTE)2_2^TEc.TEc.TT+^EE(cTT)2) 2 
2<-2 „„ 2e-3 



2< + l "t 
2 T^EE" 



TEr-TE/^EE 



(aEE(cTE)2_2a' 
2^-2 



CEE+<,TT(cEE)2-)- 



Univariate marginals 



P{Cf^\a,) 
P{CfB\ae) 



TT 

„„ 2£-3 _2«±lf> 

2?^e « 



2^-3 2^+1 '^l 



21-1 _ 21 + 1 "l^ 



Univajriate marginals, one conditional variable 



P{Cl'^\Cr,a,) 



P{Cj^\C^^,ae) 



P{Cf'^\Cj^,ae) 



P{Cr\Cf^,ae) 
P{Cf^\Cr,ae) 
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2^ — 3 
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which is a simple numerical constant. The desired 
marginal distribution therefore reads 



Note. — The determinant <t^ denotes the two-dimensional {T, E} determinant. See main text for definitions of Ii {£, A, B) and I2 {£, C, D) 
The distribution then reads 

P(Cf E|^^) ^ I p^cr,Cr\a,)dC, 

.—>..-, 9P — '^ EE 2i—l 

/ EE^2£^ _f2e+i.ze /-oo / 1 \ 



(22) 



cx 



The integral in this expression is, for f > 1, 



Li 



dy 



r(2^) ' 



(21) 



Again, wc note that the corresponding expression for 
P{Cj^\ai) is obtained simply be replacing EE with TT. 

3.4. The marginal TE distrihution, P{Cj^\ae) 

Finally, wc consider P(Cj^\(Ji) = 
P{Cj^,Cf^\ai)dCf^. As was the case for 
P{Cf'^ ,Cf^\ae), this distribution does not have a 
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closed form, but may instead be written on a computa- 
tionally convenient form, 



P(Q™|a,) cx 



TT^EE 



h{l,C,D), (23) 



where I2 denotes the integral 



2£-3 

X 2 e 



-dx. 



(24) 



The two dimensionless auxiliary parameters in this inte- 
gral are 



C 



D 




TE 



2i+l „TT „EE 

2 '^i 



(25) 
(26) 



Thus, as was the case for /i, also I2 may be tabu- 
lated over a two-dimensional grid for each multipole. It 
is therefore computationally straightforward to evaluate 
P(C™|cTf) at a given value of C™, even if it does not 
have a closed analytic expression. 

3.5. The conditional T E distribution, P{Cf^\Cf^,(Ti) 

From the above expressions, we may also derive 
all possible conditional distribution, since P{x\y) = 
P{x,y)/ P(y). Here we will only explicitly consider the 
conditional TE distribution, 



TE\r^EE 



p{ci-\c, 



p{c'r,crw,) 

P{Cf^\ae) 



cx 



2af^CP + aJTCf^)^ (^rf ^)- 



(27) 

which is relevant to several important applications (see, 
e.g.. Section [4|). 

If we make the same linear transformation of Cj^ as 
in Section [3?3l but including an additional 2£ — 2 factor. 



V2f- 2 

. /nrTriEE 



^EE^TE ^TE,-<EE\ /os\ 

- ai Cf, ) , (28) 



this we see that this may be rewritten into a familiar 
form, 

p{cr\cr,<j,) cx - — (29) 

(1+21^)^ 

We recognize this as the Student's t distribution with 
V = 2£ — 2 degrees of freedom. 

Since one of our goals is to sample from this distribu- 
tion, it is useful to have its cumulative distribution, F, 
on a closed form, 



1 ^^. ^ + l . 2Pi(i^;i; 
2+"^^ 2 > ^/^^(|; 



(30) 



Here 2F1 denotes the hypergeometric function (e.g. 
Abramowitz & Stegun 1972). 



Direct inverse Wishart sampler 
New conditional sampler 




10000 



Power spectrum, C 1(1 + \)l2iz (|j,K ) 



/ 

Fig. 1. — Comparison of marginal distributions obtained with 
the direct inverse Wishart sampler and with the new conditional 
sampler. 

4. CMB GIBBS SAMPLING APPLICATIONS 

The analytic expressions derived in Section [3] are the 
main results of this paper. Being completely general, 
these may in principle be applied to a wide range of prac- 
tical CMB applications. However, they are particularly 
useful for methods that consider both the sky signal s 
and the power spectrum Ci as free variables, such as 
the CMB Gibbs sampler. In this section, we demon- 
strate two specific Gibbs-based applications, namely 1) 
a new sampling algorithm that supports general bin- 
ning schemes, and 2) new Blackwell-Rao estimators for 
marginal and conditional distributions. 

4.1. A new Ce sampling algorithm 

The CMB Gibbs sampler draws samples from the joint 
posterior, P(s, C^jd), by alternately sampling from the 
two corresponding conditionals (e.g., Jewell et al., Wan- 
delt et al., Eriksen et al. 2004), 



■i+l 



P(s|Q,d), 
PiCi\s 



(31) 
(32) 



(In this expressions, the arrow indicates sampling from 
the distribution on the right hand side.) The former of 
these distributions is a high-dimensional Gaussian distri- 
bution, while the latter, which is the topic of the present 
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Fig. 2. — A single joint power spectrum sample Ce drawn from 
P(Ci\a(), adopting individual binning schemes for Cj"^ , Cj^ and 
^EE 

paper, is a product of independent inverse Wishart dis- 
tributions. 

There is already a well known and simple algorithm 
available in the literature to sample from the inverse 
Wishart distribution (e.g., Larson et al. 2007 or Gupta 
& Nagar 2000): Let p be the dimension of the target ma- 
trix (e.g., p = 2 for {T,E}), and = {2£ + l)at. Then 
the algorithms goes as follows: 1) Draw 2£ ~ p vectors 
from a Gaussian distribution with covariance matrix S^; 

2) compute the sum of outer products of these vectors; 

3) invert this matrix. 

This algorithm produce samples for a given multipole 
£. However, in low signal-to-noise applications it is of- 
ten desirable to bin many multipoles together, in order 
to increase the effective signal-to-noise of each variable. 
Because the CMB power spectrum essentially scales as 
0{£~^), it is customary to bin in units of C££{£ + l)/2n. 
With this convention, the above algorithm may be gen- 
eralized to include binning by redefining the covariance 
matrix as follows. 



E 



£{£ + !) 
27r 



(2£+l)(T£. 



(33) 



Here h = [^mim ^max] indicates the multipole range of the 
bin under consideration. Note that the total number of 
modes is now M = (^max + 1)^ 
Gaussian vectors must be drawn from S 



^min' ^'^'i therefore M 



Unfortunately, this method has the serious drawback 
that the binning scheme must be identical for Cj"^, Cj^ 
and Cf^. This is a problem because the signal-to-noise 
ratio of most experiments is very different for TT than 
for EE, and one would lose much spectral resolution if 
one were to bin Cj"^ with a bin size such that the signal- 
to-noise ratio for the corresponding Cf^ bin is unity. 

The new analytical expressions derived in Section [3] 
allows us to resolve this problem. First, we use the defi- 
nition of a conditional distribution, and write 



F{Lf^ ,0£ ,Cj, \ai) 



■P(Cna,). 



EE 







(34) 



The algorithm may now be written symbolically as fol- 
lows, 

(35) 
(36) 
(37) 



Cr^P{Cr\Cf^,a,) 



Then {Cj^, Cj^ ,Cf^} will be a sample drawn from the 
joint distribution P{Cj'^ ,Cj^ ,Cf^\(Ti). Note that each 
of these conditional distributions is a one-dimensional 
distribution, and the correlation with other polariza- 
tion components comes into play only conditionally, not 
jointly. One is therefore completely free to specify the 
binning schemes for each component independently of 
the others. 

We still need to write down a sampling algorithm for 
each of these conditionals. However, since these are all 
one-dimensional, this is a trivial task. The simplest ap- 
proach, and the one currently implemented in our codes, 
is to take advantage of the cumulative distribution: Sup- 
pose we want to draw a sample from a univariate dis- 
tribution, P(x), and have access to the corresponding 
cumulative distribution, F{x). Then we can draw a uni- 
form variate rj C/[0, 1], and solve for F{x) — rj. The 
sample x will then be drawn from P{x). 

However, a computational issue with this approach is 
the evaluation of the cumulative distribution. With the 
notable exception of P{Cj'^\Cf'^ ,(Ti) for a single multi- 
ple, we do not have explicit analytic expressions for the 
cumulative distributions. Therefore, in these cases we in- 
stead have to map out the analytic probability densities 
over a grid, and do the integration numerically. Fortu- 
nately, this requires only ~ 0(10^) function evaluations, 
and is therefore computationally quite fast. The cost of 
the full Gibbs sampler is anyway hugely dominated by 
sampling from the sky signal distribution P{s\Ci, d). 

Nevertheless, one might want to consider alternative 
approaches for applications in which this sampling step 
may be dominating. In such cases, the rejection sampler 
(e.g., Liu 2001) appears as a promising candidate. In 
this approach, one samples from an auxiliary distribu- 
tion that preferably should approximate the target dis- 
tribution quite well. Since our distributions are all one 
dimensional, and strongly single-peaked, it should not 
be difficult to establish such auxiliary functions. Indeed, 
the Student's t distribution itself is a typical candidate 
for such purposes, since it has a quite heavy tail. 

In Figure [1] we compare the marginal distributions ob- 
tained by the direct inverse Wishart distribution and by 



the new sampler presented here, for ^ = 10 only. As ex- 
pected, they agree perfectly. Next, in Figure[2]we show a 
single sample from the joint a\l-£ distribution P(C^|ct^), 
with appropriate binning schemes for each of Cj"^ , Cj^ 
and Cf^. Again, producing a similar sample with the di- 
rect inverse Wishart sampler is not possible, as discussed 
above. 



4.2. Marginal Blackwell-Rao estimators 

Next, we consider Blackwell-Rao estimators for 
marginal and conditional distributions, and focus for now 
on the factorization of P(C^|f7^) given in equation 1341 
However, we note that any one of the distributions in 
Table [T] gives rise to a separate estimator. 

R ecall first the der ivation of the Blackwell-Rao estima- 
tor (|Chu et al.ll2"005h . 



P{Ce\d) = I P{Ce,s\d)ds 

P{Ci\s,d)P{s\d)ds 
PiCe\ae)Piai\d) Den 

Ng 



(38) 



Thus, the full Blackwell-Rao estimator for P{Ci\d) is 
nothing but the sum (or average) of P{Ci\ae) over Gibbs 
samples, cr^. 

This estimator, as written here, has notoriously poor 
convergence properti es as the d imension of the parame- 
ter volume increases (|Chu et al.l [2005i): It requires an ex- 
ponential number of samples in order to converge. The 
reason is simply that the volume of a single Gibbs sam- 
ple is given by cosmic variance alone, whereas the volume 
of the full distribution is determined also by noise and 
sky cut. Therefore, even if the width of P{Ce\ai) for 
a single multipole is as much as, say, 90% of the width 
of P{Ce\d), the relative volume fraction in fmax dimen- 
sions is / = 0.9'™'*'^. For £max = 30, this number is 
/ = 0.04, and for ^^ax = 100 it is / = 2.65 • IQ-^. 
Clearly, a brute-force Blackwell-Rao approach will not 
work for high-dimensional spaces unless the volume ra- 
tio per multipole is unrealistically large. However, the 
same problem does not affect the marginal distributions 
described above, because the number of dimensions is 
low, and typically just one, by construction. 

Let us consider the Blackwell-Rao estimator for 
P{Cf^\d) for a single multipole, £. First, marginaliza- 
tion over other multipoles consists, as usual for MCMC 
algorithms, simply of disregarding the sample values of 
all other ts. Second, marginalization over Cj"^ and Cj^ 



Blackwell-Rao estimator 
Sampled distribution 




1000 



2000 




Power spectrum, C 1(1 + l)/2n (|iK ) 



Fig. 3. — Comparison of three Blackwell-Rao estimators and sim- 
ple histograms computed from a low-resolution CMB simulation. 



P{Cf^\d) 



is done using the expressions derived in Section [31 

//P(cr,c-,cnd,.cr.cr 
// i:/'(cr,cr,cf=i<,;)<icr<ic: 



■TE 



'TE 



r i,EE\2i^ 21+1 "l''^'^ 



(39) 



Similarly, the Blackwell-Rao estimator for 
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P{Cj^\Cf^,d) reads 

Picr\cr,d 



P(C™,C™ |d) 



(40) 



■On L//> 



(41) 



First, 



note that because Cf^ is a constant in this 
expression, P(C|'^ |d) is also a constant, and can be 
disregarded after equation 1401 Second, we emphasize 
that it is crucial to use the full joint expression for 
P(C™, C/'^jd) in this estimator, and not the naive "al- 
ternative" P(C™|cfE,d) « Y.^P^CJ^\Cf.'y\Y the 
latter approach would require the underlying Gibbs sam- 
ples, cr^, to be drawn conditionally on Cj^ in original 
Gibbs analysis, and this is usually not the case. 

Finally, we consider the expression for 
P(C7T|CTE^^EE^ J) foj. ^ gijjgig muitipole. How- 
ever, there is little simplification involved in this 
expression, as it simply reads 

p{cr\ci\cr. d) oc p[cj\ Q™, crid) 



(42) 



Y^p(cr^cr^cr\ 



E 



where at and Ct indicate 2-dimensional {T, matrices. 

As a simple illustration, we plot the three Blackwell- 
Rao estimators given above for ^ = 10 in Figure [31 as 
computed from a low-resolution simulation (A^sido = 16, 
Ciax = 47, Gaussian beam of 10° FWHM, and white 
noise of = 1/xK and a'^'^ = 0.5/^K for tempera- 
ture and polarization, respectively) drawn from a stan- 
dard ACDM spectrum. As expected, the agreement with 
direct histograms is excellent, but the smoothness and 
faster convergence of the Blackwell-Rao estimator make 
it the preferred choice for most applications. 



While these expressions are useful in their own right, 
for example for plotting marginal or joint power spec- 
tra and corresponding confidence regions from a set of 
Gibbs samples (say, by maximizing and/or integrating 
the Blackwell-Rao estimator for each £ individually), 
their main application lies in providing robust estimates 
of P{Cj^,Cj^,Cf^\d) in terms of univariates. This 
may open up for some very interesting possibilities for 
stabilizing the exponential behaviour of the full multi- 
variate estimator. This topic will be explored further in 
a future publication. 



5. CONCLUSIONS 

We have derived computationally convenient expres- 
sions for all marginals of P(Cf|(jf). These expressions 
may be useful to any CMB analysis method that consid- 
ers both the CMB sky signal s and the power spectrum 
as unknown parameters. One prominent example is 
the CMB Gibbs sampler. 

We have also presented two specific applications of 
these expressions. First, we demonstrated a new sam- 
pling algorithm for P(Ci\(T£) that supports different bin- 
ning schemes for each polarization component. This 
is useful because most experiments have very different 
signal-to-noise ratio to temperature and polarization. 

Second, we have provided explicit expressions for 
the Blackwell-Rao estimators for P{Cj'^\Cj^,Cf^,d), 
P{Cj^\Cf^,d) and P(Cf ^|d). Together, these three 
can be used to map out the joint distribution P(C£ |d) 
for a single multipole in terms of univariate distributions 
alone. Further, we note that any of the distributions 
listed in Table [1] give rise a separate, and potentially 
useful, Blackwell-Rao estimator. 
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